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1. Introduction 



Low temperature, weakly-ionized, radio-frequency 
(rf) discharges are commonly used in semiconductor 
processing. RF plasma enhanced chemical vapor depo- 
sition (PECVD) and reactive ion etching (RIE) are two 
key steps in integrated circuit manufacturing. Though 
the semiconductor industry has been recently experi- 
menting with inductive coupled plasma (ICP) and elec- 
tron cyclotron resonance (ECR) sources for high rate 
and large area processing, rf plasma reactors continue to 
be used in large-scale manufacturing. Regardless of the 
nature of the power source, all discharges used in pro- 
cessing are notoriously complex and characterized by 
several interacting phenomena: plasma generation of re- 
active species, power deposition, sheath behavior and 
ion acceleration, plasma-surface interaction, gas flow, 
and heat transfer. These phenomena seriously affect 
equipment and process design efforts. Indeed, it is a 
well-known fact that seemingly identical commercial 
reactors behave differently under "identical" operating 
conditions. This long-standing problem with plasma 
processing reactors has prompted researchers, under the 
auspices of the Gaseous Electronics Conference (GEC), 
to design a standard reactor as an experimental platform 



[1]. This GEC reactor has been used in many laborato- 
ries in the past 5 years and experimental data from 
different research groups are routinely compared to en- 
hance our understanding of rf plasma processing. The 
preceding articles in this Special Issue describe results 
from experimental studies and comparison of data taken 
from different GEC reactors. Another aim of the GEC 
reference reactor related research is to provide diagnos- 
tics data for validation of discharge models. Modeling 
of plasma processes has been intensively pursued in the 
past decade to aid in the interpretation of experiments 
and provide an understanding of discharge mechanisms. 
This article reviews one-dimensional modeling studies 
pertinent to the GEC reference reactor in the published 
literature. The number of articles on the modeling of 
GEC reactors is rather limited and this review is based 
on the following studies: Sommerer and Kushner [2], 
Lymberopoulos and Economou [3], Meyyappan and 
Govindan [4], Young and Wu [5], and Riley et al. [6]. In 
addition, the discussion is augmented with unpublished 
results from our work. For a detailed general review of 
the discharge physics and chemistry modeling, the 
reader is referred to Ref. [7]. 
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2. A Review of Models 

Direct solution of the Boltzmann equation for rf pro- 
cessing problems is difficult and not commonly done, 
even though the Boltzmann equation is the basis for all 
theoretical descriptions of the plasma. Numerical mod- 
eling of discharges falls into one of two major cate- 
gories: fluid and kinetic models. Fluid models consist of 
either two or three moments of the Boltzmann equation, 
along with Poisson's equation and provide information 
on the number density, momentum, and energy of each 
charged species. Fluid models do not solve for the elec- 
tron energy distribution function (EEDF); rather, they 
require EEDF and rate constants as input to the model. 
In contrast, kinetic models self-consistently solve for the 
EEDF and provide all of the plasma characteristics. In 
the following two sections, we provide a brief review of 
the two schemes as practiced in the studies referenced 
here [2-6]. 

2.1 Kinetic Schemes 

The most common kinetic scheme involves integra- 
tion of the Boltzmann equation with the aid of Monte- 
Carlo techniques to account for collisions. For a review 
of particle-in-cell (PIC) techniques, the reader is re- 
ferred to Ref. [8] and references therein. In PIC and 
other statistical schemes, a large number of superparti- 
cles are needed to ensure meaningful statistics and ob- 
tain smooth results. For electronegative discharges, 
where ion recombination reactions are significant, even 
a larger number of superparticles would be required. 
Hence, the kinetic schemes are computationally inten- 
sive. In contrast, fluid models require much less com- 
puter time and are more robust. Kushner and coworkers 
[2,9] introduced a hybrid scheme which retains the ad- 
vantages of both the fluid and PIC approaches and 
reduces the computational burden. Their approach is 
briefly discussed below and their results for the GEC 
reference reactor are discussed in Sec. 3.1. 

The hybrid scheme of Kushner and coworkers [2,9] 
consists of three modules, namely, an electron Monte- 
Carlo simulator (EMCS), neutral chemistry and trans- 
port model (NCTM), and self-consistent fluid model 
(SCFM). The computer simulation begins with an initial 
guess for electron density and electric field as a function 
of position and phase. The EMCS is run first using 
100-500 superparticles for 10-200 rf cycles. In EMCS, 
the trajectories of individual electron particles are 
tracked in time with prescribed time steps. When the 
number of superparticles changes due to inelastic pro- 
cesses such as ionization and attachment, simulated 
electrons are removed or added to maintain approxi- 
mately the same number of simulated particles. The 



EED is formed by periodically recording the phase 
space location of each electron in the ensemble. The 
output of EMCS is the EEDF, electron impact rate coef- 
ficients, source terms for electrons and ions, and trans- 
port data. NCTM is then used to obtain steady-state 
densities of neutrals, and source and sink rates for 
charged species. The NCTM consists of continuity 
equations for each neutral species expressed in terms of 
a finite difference formulation. The information from 
both EMCS and NCTM are next used in the self-consis- 
tent fluid model. The SCFM consists of continuity 
equations for each charged species and Poisson's equa- 
tion. The continuity equation itself comprises of drift 
and diffusion mechanisms for transport and source and 
sink terms. The equations expressed in a finite differ- 
ence form are solved using an Euler method to yield the 
charged species density and electric field at every call to 
the SCFM. The SCFM is run for 5-20 rf cycles at every 
call. The overall process between the three modules is 
iterated 100-1000 rf cycles as warranted by the time- 
scales of the problem until a periodic steady solution is 
reached. 

Riley et al. [6] also used a modified hybrid model in 
which the kinetic portion of the scheme involves direct 
solution of the Boltzmann equation for the electrons in 
one configuration-space and two velocity-space dimen- 
sions. A method of characteristics combined with a 
finite-difference expression for the collision integrals 
was used to solve the Boltzmann equation. The solution 
for the electrons alternated between the kinetic equa- 
tions and a fluid or fluid-like "average" description. The 
neutral and heavy charged particles were treated with 
fluid equations at all times. The solution was run until 
convergence. Details are contained in [6]. 

2.2 Fluid Models 

2.2.1 Governing Equations A detailed discus- 
sion on the derivation of moment equations, closure 
relations, boundary conditions and other related subjects 
can be found in Ref. [7]. Here, only a brief summary is 
given to indicate the general formulation as employed in 
Refs. [3-5]. The fluid model consists of mass, momen- 
tum, and energy conservation equations for each 
charged species j . 



dn 



! +V'¥rI^ 



(1) 



dt 



(njfrijVj) + V * (njmjVjVj) = 



-VPj + qjfijE-rijmjVjVj. 



(2) 
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dt 



Wj + V * VjWj = qjnjVj-E-V-PjVj 



J = nv = sign(q)jjuiE -DVn 



(5) 



+ V-K ] VT j -J J R l] H l . 



(3) 



In the above 



Here n is the number density, m is particle mass, and v 
is average or directed velocity. The total energy com- 
prises of kinetic energy associated with directed motion 

and thermal energy; w = n{-mv 2 +U 

equations, the sign of qj must be used: + for positive ions 
and - for electrons and negative ions. Other notations are 
as follows: t is time, R^ is the rate of generation or 
consumption of species j through inelastic collision i\ P 
is species partial pressure defined as nkT where k is the 
Boltzmann constant and T is temperature; E is electric 
field, v is elastic collision frequency; K is thermal con- 
ductivity; and Hi is energy loss following inelastic colli- 
sion process i. In addition to Eqs. (l)-(3), we need 
Pois son's equation: 



V 2 tA=- 



R 

e " 



(4) 



Here p = Xj qffi is the net charge density and 
E = -Vijj where \p is potential. 

The terms in the momentum equation [Eq. (2)] re- 
spectively represent time rate of change of momentum, 
convective acceleration or inertia, pressure gradient, 
particle drift due to electric field and finally, momentum 
loss (or drag) due to electron/gas, ion/gas elastic colli- 
sions. The momentum loss or friction term in Eq. (2) is 
written with the aid of an effective momentum transfer 
frequency v m . Equivalently, an effective momentum re- 
laxation time T m may be defined and the friction term is 
written as nmv/r m [5]. The terms in Eq. (3) respectively 
represent time rate of change of total energy, convective 
transport of energy, energy gain from the electric field, 
rate of work done by pressure forces, energy transport 
by conduction, and energy loss due to inelastic colli- 
sions. The energy loss term in Eq. (3) includes loss from 
ionization, excitation, and other inelastic collisions. Al- 
ternatively an energy relaxation term may be defined as 
in Ref. [5] and the energy loss term then is written as 
n (w-Wo)/r e where w is the energy of the background 
gas. The heavy ions exchange energy with the back- 
ground gas efficiently and their temperature is close to 
that of the gas. For this reason, rf discharge models 
commonly do not include an energy equation for ions. 

The studies in Refs. [4] and [5] use the above three 
moment equations. In contrast, the study in Ref. [3] 
approximates the momentum equation by neglecting the 
time derivative and inertia terms in Eq. (2) and uses a 
drift-diffusion relation: 



which is substituted in the continuity Eq. (1). Consis- 
tently, the kinetic energy l/2mv 2 is taken to be small 
compared to thermal energy U in Eq. (3). In Eq. (5), /jl 
is the mobility and D is the diffusivity. The set of Eqs. 
(1) and (3) with flux relation Eq. (5) form the two 
moment approach, which has one less partial differen- 
tial equation (per species) in each direction than the 
three moment approach and therefore results in consid- 
erable savings in computer time. While a two moment 
approach may be suitable to describe electron transport 
for reasonably high pressures, inertial effects are signif- 
icant in ion transport. Recognizing this, Lymberopo_ulos 
and Economou [3] used an effective electric field E in 
the flux relation [Eq. (5)] for ions where Es given by 



dt iMm \ I 



(6) 



Here subscript i denotes ions. Equation (6) has been 
derived assuming that ion transport can be modeled as 
an inertially accelerating particle: 



dVi _ F_ _qE qVj 
dt m { m x jjb\\n x 



(7) 



The net velocity v* is given by \i\E, and substituting this 
in Eq. (7) provides Eq. (6). Thus, the use of Eq. (6) 
attempts to consider ion inertia while preserving the 
simple form of Eq. (5). 

2.2.2 Rate Constants and Transport Data The 
fluid model equations presented in the previous section 
need to be augmented with information on rate expres- 
sions for various species generation/loss processes and 
transport parameters. In the kinetic schemes described 
in Sec. 1, since the EEDF is a model output, information 
on rate processes and transport parameters were gener- 
ated self-consistently using available cross-section data 
as input. This is not the case with fluid models. The 
form of EEDF has to be assumed and then the rate 
constant can be generated using known cross -sections 
from: 



Jo 



/(e)o-j(e)M(e)de . 



(8) 



Here kj is the rate constant for process j, a is collision 
cross-section, e is electron energy and u is electron 
velocity, € = 1/2 mu 2 . It is easy to see that the utility of 
any fluid model critically depends on the input for rate 
processes, i.e., the form of EEDF assumed. The 
common sources for such input are dc Monte-Carlo 
simulations, zero-dimensional (spatially homogeneous) 



465 



Volume 100, Number 4, July-August 1995 

Journal of Research of the National Institute of Standards and Technology 



solution to the Boltzmann equation, and swarm experi- 
ments. "Young and Wu [5] used a dc Monte-Carlo proce- 
dure which can be thought of as a "numerical swarm 
experiment." This involves following a swarm of test 
electrons in a specified constant electric field with colli- 
sional processes modeled statistically in a manner simi- 
lar to the kinetic scheme in Sec. 2.1. Young and Wu 
generated inelastic rate constants, mobility and diffusiv- 
ity as a function of electric field-to-gas density ratio 
(E/N). Alternatively Lymberopoulos and Economou [3] 
solved the time- and space-independent Boltzmann 
equation for a given E/N to obtain the EEDF and com- 
puted rate constants using Eq. (8). Meyyappan and 
Govindan [4] used the functional relations of rate con- 
stants vs mean energy from Ref. [3]. 

2.2.3 Boundary Conditions To complete the 
fluid model description, we need statements of 
boundary conditions for solution variables. In one- 
dimensional rf discharge models, the conditions 
discussed below are applied to both electrodes. For elec- 
trons, the net flux at the electrode is given by the sum 
of two contributions: electrons lost due to recombina- 
tion and electrons generated by secondary electron 
emission. It is usually assumed that all electrons striking 
the electrode are absorbed, i.e., unity sticky coefficient 
or zero reflection coefficient. The flux of electrons to 
the electrode is given by nv t /4 where v t is the thermal 
velocity given by ($kT e /7rm) m . Sometimes this flux is 
written k r n where k is a recombination coefficient and 
comparison with the previous expression gives k T =vJA 
[3]. Including secondary electron emission, the 
boundary condition for electrons becomes 



nvJJ 



77+ ** 



where 



q e =j e U-K 



dT 

dx ' 



(id 

(12) 



When y is zero, the electron temperature gradient is 
zero at the electrode. A constant electron temperature of 
0.5 eV at the electrode has been assumed in all of the 
previous works reviewed here [3-5]. For the three mo- 
ment approach, additional conditions are required in the 
numerical solution procedure to conveniently close the 
discretized set of governing equations. For this purpose, 
second order extrapolation conditions (zero second 
derivatives) for the electron and ion velocities are nor- 
mally used. Finally, for Poisson's equation, a zero- 
voltage condition is imposed on the grounded electrode. 
A sinusoidal voltage waveform is normally used at the 
powered electrode. 

3. Modeling Results 

It is noted that the present review is confined to one- 
dimensional simulations of the GEC cell and hence the 
only relevant dimension for modeling is the electrode 
gap, which is 2.54 cm. None of the one-dimensional 
models to date [2-6] account for the asymmetry of the 
discharge. Nevertheless, these studies have provided 
valuable information on the rf discharge mechanisms. In 
our present work, we have attempted to account for the 
asymmetry using an "effective area" approach; these 
results will be discussed in Sec. 3.2. 



7e = " 



n e v t 



-yj + 



(9) 



Here y is secondary electron emission coefficient and 
j=nv is the flux. In the two moment approach, flux 
relation [Eq. (5)] is used fory. Alternatively, if sec- 
ondary electron emission is negligible and for a unity 
sticky coefficient, electron density at the electrode is 
assumed to be zero [4, 5]. For positive ions, the flux at 
the electrode is dominantly due to drift and therefore, 



dn + /dx = 0. 



(10) 



The negative ions in rf discharges are repelled from the 
electrodes at 13.56 MHz and therefore negative ion den- 
sity at the wall may be set to zero. The boundary condi- 
tion for electron mean energy may take the form of an 
energy balance at the electrode. If the secondary elec- 
trons are assumed to emerge with an energy e s , then the 
energy balance is written as 



3.1 Hybrid Model Results 

Kushner and coworkers [2,9] modeled a number of 
gas discharges corresponding to GEC cell conditions: 
He, N 2 , He/N 2 /0 2 mixture, He/Cl 2 , He/HCl, He/CFVO, 
mixture and SiR/NHs mixture. Unfortunately, compari- 
son with measurements was limited since experimental 
results for these discharges were not available at the time 
of their work. They provided a comparison of their 
modeling results with microwave interferometry mea- 
surements of line integrated density in He from two 
GEC reference reactors. This comparison is reproduced 
in Fig. 1. The simulation results in Fig. 1 are for two 
cases, with and without inclusion of the metastables in 
the model. The electron density is seen to increase by a 
factor of 2 to 3 when contributions from Penning reac- 
tions and multistep ionization processes are included in 
the model. To compare against the microwave interfer- 
ometry measurements, the plasma was assumed to have 
no radial variation and the one-dimensional results were 
integrated over the 10.16 cm diameter of the electrode. 
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This also assumes that the plasma does not spread 
beyond the cylinder defined by the 10.16 cm electrodes. 
Despite these assumptions, the results from the 1-D 
hybrid model are in reasonable agreement with mea- 
surements, as shown in Fig. 1. The authors believe that 
at high voltages, corresponding to a power density of 
above 10 mW/cm 2 , secondary electron emission may 
become important and must be included in the model in 
order to reproduce experimental observations. Time and 
spatially averaged electron density and positive ion den- 
sity from a hybrid model of a 0.9/0.1 He/Cl 2 discharge 
are shown in Fig. 2. The plasma density exhibits a linear 
dependence on rf voltage. C1J is the dominant positive 
ion even though the mole fraction is 0. 1 . The production 
of He + is small and in addition, charge exchange with 
Cl 2 depletes the He + ions rapidly. 




100 200 300 

DISCHARGE VOLTAGE AMPLITUDE (V) 



turn 

n 

SEi 

UJ 2 

O W 

§£ 

< o 

2 tu 



10 



,12 



10 



11 



10 



10 



1#r 



~~r 



_ r 



+ MODEL WITHOUT He 
x MODFL WITH He* 
n SANDIA REFERENCE CELL 
o WPAFB REF. CELL 
1.0 Torr He 




■ 10 12 



10 



13 



10 



11 



10 



10 



tU £ 

UJ t 
> CO 

*? 2 

UJ w 

o ° 
< 2: 

CL O 
CO CE 

it 
II 



50 100 150 200 250 

DISCHARGE VOLTAGE AMPLITUDE (V) 

Fig. 1. Comparison of hybrid model results and microwave interfer- 
ometry measurements of electron density in helium discharges. 
Reproduced with permission from Ref. [2]. 



Riley et al. compared their simulation results for he- 
lium discharges with microwave interferometry mea- 
surements reported in Ref. [10] on a symmetrically 
("push-pull") driven GEC cell. This comparison for 
electron density is reproduced in Fig. 3. The compari- 
son in general is within a factor of 2 except at y rf = 150 
V and/? = 133 Pa (1.0 Torr) where the difference seems 
to be about a factor of 3. In some cases, the reported 
electron densities also seem to be in reasonable agree- 



ment with the values from Kushner and Sommerer [2] 
even though the trends appear somewhat different. Riley 
et al. also compared He metastable densities with mea- 
surements and found that though the predicted trends 
(variation with rf voltage) are as seen in experiments, the 
densities are higher than experimental data by less than 
a factor of 2 for the triplets and about a factor of 3 for 
the singlets. The authors suggest this discrepancy is due 
to the lack of detailed treatment of the high Rydberg 
levels in their model of the He excitations. 
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Fig. 2. Electron and ion densities in a He/Cl 2 discharge as predicted 
by hybrid model. Reproduced with permission from Ref. [9]. 
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in a helium discharge. Reproduced with permission from Ref. [6]. 
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3.2 Fluid Model Results 

The studies in Refs. [3-5] demonstrated the capabil- 
ity of fluid models to capture key rf discharge phenom- 
ena. These one-dimensional models used GEC cell di- 
mensions (electrode gap of 2.54 cm) and computed 
discharge behavior as a function of pressure and applied 
rf voltage. The discharge was assumed symmetric and 
all radial variations were ignored. No comparison 
against GEC cell measurements were attempted in these 
early models; for this reason, only a brief summary 
highlighting key achievements of each work is presented 
here. Lymberopoulos and Economou [3] were the first 
to include a detailed balance equation for metastables in 
fluid simulations of rf discharges. At a pressure of 133 
Pa, their results showed that the electron density with 
metastables included in the model is about six times 
higher than that without metastables. Stepwise-ioniza- 
tion was shown to be the dominant mechanism for elec- 
tron production at 133 Pa compared to direct ionization 
of the ground state atoms. Their results emphasized that 
neutral transport and reactions must be considered in a 
self-consistent manner, even in inert gas discharges, in 
order to capture the observed behavior of rf discharges. 
Lymberopoulos and Economou also effectively demon- 
strated an approach based on the zero-dimensional solu- 
tion of the Boltzmann equation to generate the rate con- 
stants needed in fluid models. This approach is more 
realistic than the Arrhenius expressions based on 
Maxwellian EEDF commonly used in fluid models until 
that time. Young and Wu [5] showed that the three- 
moment fluid model is much more capable of capturing 
non-equilibrium effects than one- and two-moment 
models. Comparing their results against Monte-Carlo 
simulations, they concluded that the three-moment fluid 
model can predict, within reasonable accuracy, the tran- 
sition between the a and y regimes. Meyyappan and 
Govindan [4] showed that with rate constants obtained 
from a zero-dimensional solution of the Boltzmann 
equation, their three-moment model reproduces excita- 
tion wave forms as seen in experiments. They also 
showed that the three-moment fluid model can repro- 
duce several critical discharge features, previously re- 
ported only by PIC simulations, such as mechanisms of 
electron heating and relation between electron power 
deposition and ionization. The simulations in Ref. [4] 
provided some scaling laws for the discharge behavior. 
Electron density in electropositive discharges scales 
nearly linearly with rf voltage. Both the rf and ion cur- 
rent also scale linearly with rf voltage. Power deposited 
to the electrons varies linearly with rf voltage while ion 
power varies linearly with the square of the rf voltage. 

In the remainder of this section, we discuss unpub- 
lished results from our research on GEC reference reac- 



tor modeling. These results are for a 13.3 Pa (100 
mTorr) argon discharge at three zero-to-peak (V rf ) 
voltages of 75 V, 100 V and 150 V. The results are 
compared with the current-voltage, microwave interfer- 
ometry, and optical emission measurements reported in 
Refs. [10-13]. It is noted that Greenberg and Hebner 
[10] powered both electrodes in a push-pull fashion and 
their discharge was fairly symmetric while other works 
in Refs. [11-13] reported an asymmetric discharge. 
Since 1-D models are naturally symmetric, we first 
compare the results of our 1 -D model with the results of 
Greenberg and Hebner [10] in Figs. 4 and 5. Figure 4 
shows the zero-to-peak (7 rf ) current at 13.3 Pa from our 
simulations and Ref. [10]. An area of 81 cm 2 was used 
to convert the peak current density to peak current 
[6,10]. The comparison between simulation and mea- 
surement is reasonable and well within a factor of two. 
Figure 5 shows the peak electron density from the sim- 
ulation and microwave interferometry results reported 
by Greenberg and Hebner. The line averaged density 
was converted in Ref. [10] by assuming a uniform elec- 
tron density confined within the 10 cm of the elec- 
trodes. Here, the model results are nearly eight times 
smaller than the measurements. We notice that at low 
pressures, there is a significant difference between mi- 
crowave interferometry and Langmuir probe results. For 
example, such a comparison by Overzet and Hopkins 
[11] shows that at 13.3 Pa and an applied peak-to-peak 
voltage of 230 V, the interferometry result is higher than 
Langmuir probe density by a factor of 2.5. It is possible 
that the plasma may not be confined within the elec- 
trodes at 13.3 Pa but this alone cannot explain the devi- 
ation. The model does not account for secondary elec- 
tron emission, which may explain the underprediction of 
density to some extent. In any case, further simulation 
work is needed for a range of pressures and voltages. 

Next we discuss 1-D modeling results for an asym- 
metric discharge. The governing equations for a 1-D 
discharge have been extended to include variation of the 
area of the discharge in asymmetric discharges due to 
unequal electrode areas. In the case of the GEC cell, the 
walls of the cell are grounded and effectively increase 
the area of the grounded electrode. This effect was in- 
cluded in the present study by (a) increasing the area of 
the grounded electrode and by (b) assuming the area of 
the discharge varies linearly from the area of the pow- 
ered electrode to the effective area of the grounded 
electrode. The area of grounded surfaces in the GEC 
cell is approximately 34 times the area of the powered 
electrode. From available experimental data, clearly this 
is too large an area to be included as is in the effective 
area of the grounded electrode. All elements of the 
grounded surface do not contribute equally to increasing 
the effective area of the grounded mode electrode. 
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In the present study the area of elements of the grounded 
surface were weighted inversely by the minimum dis- 
tance from the element to the powered electrode. This 
choice of weights for the elements of the grounded sur- 
face was based on the assumption that the influence of 
the elemental area of the grounded surface is directly 
proportional to the mean electric field between the ele- 
ment and the powered electrode and this mean electric 
field varies inversely with the distance from the powered 
electrode. Thus, 



'-JG 



dA 



where d is the electrode gap, / is the minimum distance 
between the elemental area and the powered electrode 
and the integration is carried out over all grounded sur- 
faces. An effective area of approximately six times the 
electrode area was obtained for the GEC cell. Areas of 
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Fig. 4. Current-voltage characteristics at 13.3 Pa predicted by a 
three-moment fluid model. Measurements from Ref. [10] for a sym- 
metric discharge are also shown. 
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Fig. 5. Peak electron density from a three-moment fluid model at 
13.3 Pa. Microwave interferometry measurements from Ref. [10] for 
a symmetric discharge are also shown. 



ports in the cell were excluded from the estimate of the 
effective area of the grounded electrode. 

Figure 6 shows the phase profile of electron density 
at 13.3 Pa and V rf = 100 V. We notice that the peak 
occurs off-center toward the powered electrode, as seen 
in Ref. [12]. In an asymmetric discharge, the densities 
near the grounded electrode are low throughout the 
cycle. There is a strong modulation of density near the 
powered electrode. The peak value predicted here is 
nearly an order of magnitude different from the mea- 
surements in Refs. [11,12], as discussed earlier. Our 
model has been bench-marked against simulation results 
from kinetic schemes and measurements from Godyak 
for a helium discharge in a workshop organized at the 
1993 Gaseous Electronics Conference [13]. Further 
work at various pressures and bias levels is needed to 
draw conclusions. Figure 7 shows the corresponding 
electron temperature profile which is also asymmetric. 
The bulk of the discharge exhibits an electron tempera- 
ture of ~4 eV. The metastable (3/? + 3p0 density is 
plotted in Fig. 8. The peak value is approximately 
3X10 17 m" 3 . The metastable density does not show sig- 
nificant modulation in time. Finally, optical emission 
due to the 750.4 nm argon line is shown in Fig. 9. The 
experimental data from Ref. [14] is also shown for com- 
parison. In modeling this process, a number density 
balance for the 2/?i excited state was included. Only 
excitation from the ground state was included using rate 
constants given by Ferreira et al. [15]; production of 2/?i 
excited states through intermultiplet transitions is 
smaller by an order of magnitude. The Einstein coeffi- 
cient for the emission process (2/?i — » ls 2 ) was taken to 
be 0.445 X10 8 s _1 . We note from Fig. 9 that the model 
predicts the spatial location of the peak very well, indi- 
cating that our effective area approach is reasonable. 
Conventional, symmetric 1-D models would have pre- 
dicted two spatial peaks displaced symmetrically about 
the center of the discharge, one at t = 0.25 and another 
at t = 0.75 [3,4]. The model prediction of emission 
intensity is within a factor of 2-3 of the experimental 
data. The only deviation is a band of somewhat stronger 
emission at around t = 0.25 when the voltage is at its 
peak. In the experiment, a dc bias is present across the 
electrodes, which would tend to decrease the emission at 
t = 0.25. This bias was not included in the simulations. 

4. Summary 

We have presented a review of one-dimensional 
models of the GEC reference plasma reactor based on a 
small number of publications in the literature. Most of 
the studies reviewed were limited to qualitative compari- 
son with available data and the models reproduce many 
of the observed features. Quantitative comparison be- 
tween model predictions and measurements has been 
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Fig. 6. Phase plot of electron density in an argon discharge from a three-moment fluid model, p = 13.3 
Pa, V rf =100 V. 
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Fig. 7. Phase plot of electron temperature in an argon discharge from a three-moment fluid model 
p=13.3 Pa, Vrf=100 V. 
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Fig. 8. Phase plot of metastable density in an argon discharge from a three-moment fluid model, p = 13.3 Pa, V rf = 100 
V. 
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done in a few cases in terms of current-voltage and 
line-averaged electron density; here the agreement 
ranges from reasonable to within an order of magnitude. 
The 1-D model has also been modified to account for 
the asymmetry of the GEC cell using an effective area 
approach which seems to predict the spatial asymmetry 
of emission profiles reasonably. Modeling techniques, 
both fluid and kinetics schemes, have improved vastly in 
the past 5 years. Based on the number of papers at the 
1994 Gaseous Electronics Conference reporting on di- 
agnostics in the GEC cell and the preceding articles in 
this Special Issue, it appears that much more experimen- 
tal data for a number of gases will be available for model 
validation. This provides an ideal opportunity for the 
modelers to fine tune their discharge physics modules 
and provide a basis for integration with reactor models. 
There are several challenges ahead in GEC cell model- 
ing: modeling electronegative discharges with proper 
accounting of various inelastic processes, coupling 
plasma transport equations with external circuit 
boundary conditions, investigation of the sensitivity of 
model results to the boundary conditions, and efficient 
two-dimensional time-dependent models that would al- 
low routine simulation of GEC and commercial reactors. 

Note added in proof 

The predicted electron density values (Figs. 5 and 6) 
from the three-moment fluid model increase by a factor 
of two if the ion mobility variation with the electric field 
is included in the model. The model results then are 
within a factor of two of the microwave interferometry 
measurements at pressures of 33.3 Pa and above. At 
13.3 Pa, the model predictions now are smaller than the 
experimental results [10] by a factor of four. However, 
microwave interferometry measurements at 13.3 Pa in 
argon from various groups exhibit a considerable 
disagreement, as reported by Overzet in this Special 
Issue. 
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